Vacancy effects in an easy-plane Heisenberg model: reduction of T c and 

doubly-charged vortices 
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Magnetic vortices in thermal equilibrium in two-dimensional magnets are studied here under the 
presence of a low concentration of nonmagnetic impurities (spin vacancies). A nearest-neighbor 
Heisenberg (XXZ) spin model with easy-plane exchange anisotropy is used to determine static ther- 
modynamic properties and vortex densities via cluster/over-relaxation Monte Carlo. Especially at 
low temperature, a large fraction of the thermally generated vortices nucleate centered on vacancies, 
where they have a lower energy of formation. These facts are responsible for the reduction of the 
vortex-unbinding transition temperature with increasing vacancy concentration, similar to that seen 
in the planar rotator model. Spin vacancies also present the possibility of a new effect, namely, the 
appearance of vortices with double topological charges (±47r change in in-plane spin angle), stable 
only when centered on vacancies. 

PACS numbers: 75.10Hk, 75.30Ds, 75.40Gb, 75.40Mg 



I. INTRODUCTION 

The vortex-unbinding transition in two-dimensional 
(2D) spin models with planax. symmetry (Berezinskii- 
Kosterlitz-Thouless transitionldum) has attracted interest 
recently with respect to the influence of nonmagnetic im- 
purities or spin vacancies in the lattice. In any real phys- 
ical system, some fraction of the atoms could be substi- 
tuted by impurities, and if these are nonmagnetic, the 
spins neighboring the impurity will be strongly affected 
by the missing exchange interactions. Not only could 
missing bonds cause locally lower energy densities, but 
they give the neighboring spins more freedom of motion, 
which can be expected to increase the local spin fluctu- 
ations. This can be expected to affect the static con- 
figurations, the thermal equilibrium properties, and the 
dynamic correlations, such as in EPR measurements on 
antiferromagnets 

Significant vacancy effects on the static vortex (or an- 
tivortex) configurations of ferromagnets (and antiferro- 
magnets with two sublattices) have already been found 
for a 2D easy-plane Heisenberg model (three spin compo- 
nents). Zaspel et al.U found that the critical anisotropy 
strength (5 C = 1 — A c , see Hamiltonian below) needed to 
stabilize a vortex in the planar configuration on a square 
lattice is reduced from S c i=a 0.2966 to the much lower 
value S cv i=a 0,0429 when the vortex is centered on a va- 
cancy. WysinQ found a similar result at higher precision 
(5 CV i=a 0.0455), and determined that a vacancy at the 
center of a circular system with free boundaries produces 
an attractive potential for a vortex. Using dynamic relax- 
ation and Monte Carlo simulations, Pereira et alB found 
that a single vacancy in a square system with antiperi- 
odic bounary conditions provides an attractive potential 
for a vortex. These works demonstrated a significant 
energy reduction for a vortex formed on a vacancy, com- 
pared to one formed in the center of a cell of the lattice, 
whose value depends on the type of lattice and the bound- 
ary conditions. The resulting vortex-on-vacancy binding 



energy was found to increase with increasing easy-plane 
anisotropy strength. Roth analytic and numerical calcu- 
lations by Paula et aln show that holes cut out of a spin 
lattice similarly produce interesting attractive effects on 
vortices. 

Continuum model . calcu lations for the closely related 
planar rotator modeKMiil were interpreted to suggest a 
repulsive potential between a planar vortex and a non- 
magnetic impurity however, this seems contradictory to 
later calculations. Studies of^ 2D isotropic Heisenbeca 
antifcrromagnet by M61 et al£3 and Pereira and Piresll^l 
found oscillatory dynamic modes of solitons pinned to va- 
cancies, confirming the presence of an attractive restoring 
potential. Conadering-these most recent calculations for 
several models,0ffl'E3Ei3 in general it has been seen that a 
spin vacancy attracts vortices (or antivortices) and lowers 
their energy of formation. 

In terms of the equilibrium thermodynamics, the ef- 
fect of a concentration of vacancies on the BKT transi- 
tion temperature T c of the easy-plane Heisenberg model 
has not been studied. On the other hand, Leonel et alM 
performed Monte Carlo (MC) simulations of the planar 
rotator model (two-component spins) and found a low- 
ering of the transition temperature with increasing va- 
cancy density. It was argued that vacancies produce an 
effective repulsive potential for vortices, thereby increas- 
ing the nucleation of pairs and lowering the transition 
temperature, but the vortex density was not measured 
in the MC simulations. Using the helicity modulus to 
determine T c , they found that T c goes to zero when the 
vacancy concentration of a square lattice reaches about 
30%. A similar lowering q£-T c also appears in the MC 
simulations of Berche et aZo for the same model, deter- 
mined by fitting the exponent of the spin-spin correlation 
function to the critical point value, r\ — 1/4. These lat- 
ter authors found that T c did not fall to zero until the 
vacancy concentration reached 41%, a number related to 
the percolation threshold for a square lattice. In-Sr-related 
bond-diluted planar rotator model, Castro et a/.E-3 used a 
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self-consistent harmonic approximation with vortex cor- 
rections, determining the reduction of T c with dilution, 
and the temperature variation of the correlation function 
and its exponent 77. 



Here it is interesting to consider whether a similar 
vacancy-induced reduction of T c occurs in the anisotropic 
Heisenberg model, which is a more realistic spin model 
that has a true time dynamics. An analysis of the vor- 
tex densities in thermal equilibrium, in the presence of 
vacancies, helps to explain the role of vacancies in gen- 
erating spin disorder around the transition temperature. 
The vortices in the anisotropic Heisenberg model also can 
be expected to have planar or out-of-ipj ane gtsucta m. de- 
pending on the anisotropy strength. llaoMyE°E2lG3 At 
stronger anisotropy [i.e., for the XY model, A = 0, see 
Eq. (Q)], the stable static vortices are planar, whether 
pinned on vacancies or free from the vacancies B0 Al- 
terntively, at weak anisotropy, both the stable pinned 
and free vortices have nonzero out-of-plane spin compo- 
nents, which might be expected to significantly modify 
some equilibrium properties as well as dynamic correla- 
tions. Therefore, here we present MC simulations for 
three different anisotropies, calculating the changes in T c 
and the behavior of the vortex densities, as well as other 
thermodynamic properties. 



It has been customary only to search for singly charged 
vortices appearing in MC simulations of pure easy-plane 
spin systems. Looking in individual unit cells (plaque- 
ttes) of the lattice, a net rotation of the in-plane spin 
angles through ±2tt as one moves around the cell indi- 
cates the presence of a singly-charged vortex (q = ±1). 
When vacancies are present, however, the searching for 
vorticity must be modified. Here, we searched for net 
vorticity also in the four unit cells surrounding any va- 
cancy of the square lattice. This allows for the appear- 
ance of a new effect, namely, the presence of q = ±2 
vortices, which always form centered on the vacancies. 
They appear as a very small fraction of the total vortic- 
ity density, and are present regardless of the anisotropy 
strength. Apparently, by pinning on vacancies, q = ±2 
vorticies lower their energy sufficiently due to the missing 
spin site, leading to greater ease in their thermal forma- 
tion. In addition, at low temperatures, it is found that 
most vortices (either q = ±1 or q = ±2) form initially 
on the vacancies, which gives an interesting view of how 
vacancies modify and even control the BKT transition. 



After further definition of the model, we describe the 
MC simulations, determinations of T c using finite-size 
scaling of the in-plane susceptibility, and the vacancy ef- 
fects at various anisotropies. This is followed by some 
preliminary analysis of the stability properties of the 
doubly-charged vorticies. 



II. EASY-PLANE MODEL WITH RANDOM 
REPULSIVE VACANCIES 

The model to be investigated has classical three- 
component spins defined at the sites n of a 2D square 
lattice with unit lattice constant. The spins can be ana- 
lyzed either in terms of their Cartesian components or us- 
ing polar spherical coordinate angles, S = (S x 7 S y 7 S z ) — 
S (sin 6* cos sin sin 0, cos 6*). The system is an L x L 
square with periodic boundary conditions. We consid- 
ered L ranging from 16 to 128, using the dependence of 
the thermodynamic averages on L to get estimates of the 
critical temperature in the infinite size system. 

A small vacancy density p mc is introduced into the lat- 
tice as follows. An occupation number p n for each site 
is set to the static values 1 or depending on whether 
the site n is occupied by a spin or is vacant. The frac- 
tion /3 vac of the sites has p n set to zero. (Equivalcntly, 
one can keep the spins at the vacant sites but set their 
lengths to zero.) In order to have the most simplified sit- 
uation, the vacant sites are chosen randomly, but no two 
are allowed to be within the second nearest neighbor dis- 
tance of -\/2 (the diagonal separation across a unit cell of 
the lattice). In this way, the immediate neighborhoods 
of all vacancies are equivalent: each vacant site is sur- 
rounded by eight occupied sites. This condition greatly 
simplifies the algorithm for searching for localized vor- 
ticity around the vacant sites. On the other hand, it 
limits the possible density of vacancies to be less than 
0.25 of the lattice sites (achieved in the ordered configu- 
ration having alternating rows of the lattice fully and half 
occupied by spins). In actual practice, by choosing the 
vacant sites randomly and enforcing this constraint (i.e., 
quenched repulsive vacancies), the maximum achievable 
vacancy density is p V ac ~ 0.1872 . As a result, a vacancy 
density needed to push the BKT transition temperature 
down to zero cannot be achieved, and we do not consider 
this aspect of the model here. Instead, we are more inter- 
ested in the role the vacancies play in controlling where 
the vortices are forming. 

Nearest neighbor unit length spins (S — 1) in this 
model interact ferromagnetically (exchange constant J > 
0) according to a Hamiltonian with easy-plane anisotropy 
specified by parameter A, 
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(1) 

The XY- model results for A = 0. Values of A below 1 
describe a system where z is the hard axis and xy is the 
easy plane, allowing for the appearance of vortices. The 
total number of spins in the system is 



N = N n 



(1 



)L 2 



(2) 



In general, calculated thermodynamic quantities are 
quoted here as per-occupied-site average values, i.e., nor- 
malized by N . 
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III. MC SIMULATIONS 

Classical Monte Carlo algorithms were used to esti- 
mate static thermodynamic quantities as functions of 
temperature T, with emphasis on the internal energy 
e = E/N, specific heat c = C/N, and magnetic sus- 
ceptibility of the in-plane spin components, x, all per- 
occupied-site quantities, as well as the vorticity densities 
per occupied site. It is understood that a certain percent 
of vacancies p vac has been produced in the L x L lattice 
under study, at randomly selected positions as described 
above. We found there to be very little variation in the re- 
sults with the choice of equivalent systems with different 
vacancy positions, especially for the larger lattices (i.e., 
a large system is self-averaging) . Therefore, no averaging 
over different systems at a given L was performed. 



are needed to implement the Wolff algorithm. It means 
that the Wolff clusters being formed could actually span 
across vacant sites. Clearly, this means that a large clus- 
ter being formed might actually be composed from sev- 
eral sub-clusters connected by vacant sites, a situation 
that probably enhances the mixing produced by the al- 
gorithm. 

For a single MC step (an MC pass through the lat- 
tice) , we attempted N over-relaxation moves, followed by 
N single-spin moves, followed by N cluster moves. An 
initial set of 5000 MC steps was used to equilibrate the 
system. The averages shown here result from a sequence 
of 300,000 MC steps at each individual lattice size and 
temperature. For most of the data, the error bars are 
smaller than the symbols used, hence, error bars have 
not been displayed. 



A. MC Algorithm 

The MC techniques used here have been described in 
Ref. ^ and are based partly on simulation methods de- 
veloped in Refs. |2^ , |25| . |26| . p7[ We applied a combina- 
tion of-Metropolis single-spin moves and over-relaxation 
movesES that modify all three spin components, and in 
addition, Wolff single-cluster operationso'El that modify 
only the xy spin components. The single spin moves and 
over-relaxation moves were applied to sites selected ran- 
domly in the lattice; similarly, the initial sites for cluster 
generation were selected randomly. 

In the single-spin moves, randomly selected spins were 
modified by adding small increments in random direc- 
tions, and then renormalizing the spins to unit length, 
accepting or rejecting each change according to the 
Metropolis algorithm. 

The over-relaxation and cluster moves are important 
at low temperatures, where the xy spin components tend 
to freeze and single spin moves become inefficient. Over- 
relaxation and cluster moves have the tendency to change 
spin directions with no or very small changes in energy, 
hence, their efficacy at low temperature. 

The over-relaxation moves used here consist of reflect- 
ing a randomly selected spin across the effective magnetic 
field due to its neighbors, 

Bn = jY,Pn+4SZ+J + Sl+J + XSZ+J], (3) 
a 

while preserving the spin length. All spin components are 
involved in the process, and the z components become 
more greatly involved when the anisotropy parameter A 
approaches 1. This spin change exactly conserves the 
energy, while effectively mixing up the spin directions. 

The Wolff cluster algorithm (and computer subrou- 
tine) used here is identical to that used for the pure sys- 
tem without vacancies. In the actual computations, the 
spins of the vacant lattice sites are set to zero length 
(equivalent to setting occupation p n — 0), and the cal- 
culations proceed normally. No other significant changes 



B. MC Measurements 

In terms of temperature T and Boltzman's constant k, 
the system's thermodynamic energy E and heat capacity 
C are defined via usual relations, 

E=(H), C = k[(H 2 ) - (H) 2 }/T 2 . (4) 

The instantanteous total magnetization of the system is 
the sum over all spins 

M = 5>„Sn. (5) 

For purposes of finding T c , it is important to calculate the 
associated per-spin susceptibility \ aa of any component 
a, derived from the magnetization fluctuations, 

X aa = {{Ml) - (M a f)/(NT). (6) 

Both \ xx and x vy were computed by @ and then aver- 
aged to get the in-plane susceptibility, 

X=(X xx +X vv )/2- (7) 

Finite size scaling of x was found to be the best method 
to determine T c precisely, see below. 

In the thermodynamic limit, according to the Mermim- 
Wagner theorem, (M) — > at any temperature, and this 
holds in an approximate sense in the MC averages of fi- 
nite systems. Therefore it is also interesting to calculate 
the system's total in-plane absolute valued magnetic mo- 
ment (order parameter M*), which only tends to zero in 
the high-temperature phase, and its associated per-spin 
susceptibility x*, 

M* = {^M 2 +M 2 ), X * = [(M 2 + M 2 )-M* 2 ]/(NT). 

(8) 

Related per-spin energy, specific heat, and order param- 
eter (e, c, m*), are obtained by dividing each by the num- 
ber of occupied sites, N. 
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FIG. 1: For the model with edge L = 64, at the XY limit, the 
internal energy, absolute magnetization, and specific heat per 
spin for the uniform system and with 16% vacancy density. 



For example, at A = 0, L = 64, typical results for 
the energy, absolute in-plane magnetization and specific 
heat per spin are shown in Fig. [l], comparing the pure sys- 
tem with that at 16% vacancy concentration. Note that 
the energy and specific heat per spin have rather weak 
dependence on the system size L, while m* acquires a 
sharper dropoff with inceasing L. The most obvious ef- 
fect of p vac > is the lowering of the BKT transition 
temperature. A less obvious effect is the lowering of the 
per-spin energy and specific heat in the high-temperature 
phase. This quite possibly results because a large frac- 
tion of the vortices produced in the high-temperature 
disordered phase are localized on the vacancies, as found 
below. When thus formed, vortices require a lower nu- 
cleation energy, and the system can reach a specified en- 
tropy at a lower overall energy cost. 



C. Critical Temperature 

Initially, the fourth order .in-plane magnetization cu- 
mulant Ul due to BindeiBfEl was calculated to aid in 
location of the transition temperature in the thermody- 
namic limit. It is defined using a ratio, 



U L = 1 



{{Ml 



Ml?) 



mi 



M1Y 



(9) 



This quantity becomes 0.5 in the low-temperature or- 
dered limit, and tends towards zero in the disordered 
high-temperature limit. When measured at the critical 
temperature, its value is expected to be approximately 
independent of the system size. Therefore, T c can be es- 
timated by plotting Ul vs. T for different system sizes 
and observing the common crossing point of the data. 
This definition of Ul is analogous to the more familiar 
form that would be applied to a single in-plane spin com- 




FIG. 2: Application of the //-dependence of the fourth order 
cumulant Ul on various system sizes to estimate T c / J fs 0.48 
(common crossing point of the data) at 16% vacancy density 
in the XY model. 



ponent or single-component model, viz. 

{M% ) 



ui x) = 1 



3(M2 



(10) 



Example application of Ul for finding T c for the XY 
model at 16% vacancy concentration is shown in Fig. |^. 
The transition temperature is lowered to T c w 0.48J, 
considerably less than T c w 0.70J that holds at zero va- 
cancy concentration. 

It is seen, however, that Ul requires an excessive 
amount of calculations even to. get two-digit precision for 
T c . Following Cuccoli el al£3 and their analysis of the 
pure XXZ model, a finite scaling analysis of the in-plane 
susceptibility \ is seen to be much more precise and ef- 
ficient for finding T c . The essential feature needed here 
is that near and below T c , the susceptibility scales with 
a power of the system size, 



X oc L 



(11) 



where the exponent r\ describes the long distance behav- 
ior of in-plane spin correlations below T c , see Ref. ^ for 
details. Importantly, at the transition temperature for 
the XY model, one has r\ — 1/4. Here we make the as- 
sumption that 77 = 1/4 at T c also for the models with 
A > and with vacancies present. The validity of this 
assumption is partially tested by the quality of the scal- 
ing that it produces. 

Using 77 = 1/4, we plotted x/L 7 ^ 4 versus T for the data 
from different system sizes, L = 16,32,64,128 together 
on one graph. The common crossing point of the curves 
locates the critical temperature, for example, the result 
T c / J » 0.699±0.001 is easily reproduced for the vacancy- 
free XY model. An example of this is given in Fig. |^, for 
A = at 16% vacancy concentration. An exceptionally 
tight crossing point occurs at the critical temperature, 
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FIG. 3: Application of the finite-size scaling of in-plane sus- 
ceptibility to estimate T c /J m 0.478 (common crossing point 
of the data) at 16% vacancy density in the XY model, using 
exponent rj = 1/4. 



T c /J w 0.478 ± 0.001. The clarity of the crossing point 
gives considerable confidence in the r\ = 1/4 assumption, 
even when vacancies are present. Similar results hold for 
the othe r models studied (nonzero A and nonzero p v 
see 



[HE ) where the scaling estimates of T c give dramatic 



improvement upon the more approximate estimates using 
Ul, from the same MC data. 
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FIG. 4: A spin configuration from MC simulations for L — 32, 
A = 0, p V ac = 0.16, at T = 0.85 J, with vortices indicated 
by ± signs. The projections of the xy spin components are 
shown as arrows, with line and triangular heads indicating 
positive/negative z spin components. The three larger plus 
signs are vortices of charge q — +2 centered on vacancies. 
Many vortices have formed centered on vacancies. 



D. Vortex densities 

In a system with vacancies, the presence of unit 
charged and doubly charged vortices is determined as 
follows. 

If a unit cell or plaquette is found to be fully occupied 
by spins, then the vortex search takes place in the usual 
way, counting the net vorticity there by summing the in- 
plane angular changes around the cell and normalizing 
by 2tt: 



- E 



bond ■ 



(12) 



edge bonds 



It is understood that each difference between two in-plane 
spin angles along one edge segment must be taken on the 
primary branch: — ir/2 < A^bond < 7t/2. Then q within 
a cell is forced to be an integer. In practice, the possible 
outcomes for q are and ±1, as higher charged vortices 
are unstable within a single cell of the lattice, and never 
occur in Monte Carlo simulations. 

Additionally, the search for vorticity can also be per- 
formed easily around the quartet of unit cells that sur- 
rounds an individual vacancy. A vacancy is surrounded 
automatically by eight occupied sites, connected by eight 
bonds (under our assumption of repulsive vacancies). 
Then again Eq. ( |l2] ) can be applied to determine the to- 
tal vorticity within these four cells nearest the vacancy, 



summing over the in-plane angular changes in all eight 
bonds. Now it is seen that the result for q can take the 
additional possible values q = ±2, i.e., doubly charged 
vortices are found to be stable entities when localized 
on the vacancies, but never are found to occur separated 
from a vacancy. 

An example of a state with doubly charged vortices 
is given in Fig. ||, produced in the MC simulations with 
L = 32, A = 0, p V ac = 0.16, at T = 0.85J, well above 
the critical temperature (T c w 0.478J) for this vacancy 
concentration. The locations of the q = 2 vortices are in- 
dicated by the larger plus signs; two near the top-center 
and one in the lower-left section of the system. Other 
singly charged vortices are indicated by the smaller ± 
signs. One can also note the considerable number of vor- 
tices (of any charge) that form exactly centered on the 
vacancies. 

For a state in which there are n\ singly charged vortices 
(q either +1 or -1) and ri2 doubly charged vortices (q 
either +2 or -2), the total absolute vorticity density was 
defined relative to the occupied spin sites, and giving a 
double weight to the double charges, 



P = 



(13) 



N N 

Additionally, the vorticity fraction /dbi that corresponds 



6 



0.2 



0.15 



0.05 



Vortex densities 


- 


1=0, L=64 




0.1 xf . (0.16) 


J /p(0) 







0.5 



T/J 



1.5 



FIG. 5: Thermally induced vorticity density for the uniform 
XY model [p(0)] and at 16% vacancy density [p(0.16)]. Also 
displayed are the vorticity fraction pinned on vacancies [/ P m] 
and the fraction with doubled charges [/aw], both when p V ac = 
0.16. 



to doubly-charged vortices was tracked, 

2n 2 



dbl 



(14) 



Indeed, both the q = ±1 and q = ±2 vortices are 
commonly found centered on the vacancies. Therefore, 
we also calculated the fraction / p j n of the total absolute 
vorticity that is found centered on vacancies, or, pinned 
on the vacancies: 



/pin — 



(pinned) i 



(15) 



where the sum in the denominator is over all vortices 
found in the system. As already mentioned above, the 
doubly charged vortices are always found pinned on the 
vacancies. Furthermore, at low temperatures with very 
low vortex density, essentially all vortices nucleate on va- 
cancies. 

Typical results for these various vorticity densities in 
the XY model at L = 64 are shown in Fig. g. Considering 
the curves for 16% vacancy concentration, it is significant 
that for temperatures near T c , the pinned vorticity frac- 
tion is around 75%. This is reasonable, because pinned 
q = 1 vortices have considerably lower energy than free 
ones and therefore will dominate at the lower tempera- 
tures. On the other hand, doubly-charged vorticity does 
not appear with significant population until well into the 
high-temperature phase, when it composes up to several 
percent of the total vorticity in the system. 



E. Variations with A 

The previous sections presented vacancy effects in the 
XY model, A = 0. MC simulations were also carried 



TABLE I: Dependence of critical temperature T c (p vac ) on 
anisotropy constant A, for the pure model (p vac = 0) and at 
Pvac = 0.16, obtained by the scaling of in-plane susceptibility. 



A 


T c (0)/J 


T c (0.16)/J 


0.0 


0.699 ±0.001 


0.478 ± 0.001 
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FIG. 6: For the model with edge L — 64, at the vortex- 
in-plaquette critical anisotropy, the internal energy, absolute 
magnetization, and specific heat per spin for the uniform sys- 
tem and with 16% vacancy density. 



out at two nonzero values of the anisotropy parameter: 
1) the vortex-in-plaquette critical anisotropy (A c = 0.70) 
and 2) the vortex-on-vacancy critical anisotropy (A cv = 
0.9545). At A c large out-of-plane magnetization fluctua- 
tions might be expected if free vortices were dominating 
the dynamics. At A cv large out-of-plane magnetization 
fluctuations might be expected if vortices pinned on va- 
cancies were dominating the dynamics. 

At these nonzero A, the effects due to vacancies are 
similar to those found at A = 0: reduction of T c , signifi- 
cant fraction of pinned vorticity in the low-temperature 
phase, and appearence of doubly charged vorticity in the 
high-temperature phase. 

These limited results for Tr as determined by scaling 
of x are summarized in Table |l|. At 16% vacancy concen- 
tration, the general dependence of T c on A mimics that 
found for the pure model; T c changes very little until A 
becomes very close to 1. 

The per-spin energy, absolute in-plane magnetization, 
and specific heat at A = A c are shown in Fig. |6|, where a 
mildy different result is seen compared to the XY model. 

At A cv , stronger effects are found, as seen in Fig. [t]. 
The transition temperature is reduced to T c /J ps 0.404 
when 16% vacancies are present, Fig. ||, compared to 
T c /J sa 0.608 for the pure system. The vorticity den- 
sity results are shown in Fig. O, and mimic those found 




FIG. 7: For the model with edge L = 64, at the vortex- 
on-vacancy critical anisotropy, the internal energy, absolute 
magnetization and specific heat per spin for the uniform sys- 
tem and with 16% vacancy density. 



FIG. 9: Thermally induced vorticity density [p(0)] at the 
vortex-on-vacancy critical anisotropy with 16% vacancy den- 
sity [p(0.16)]. Also displayed are the vorticity fraction pinned 
on vacancies [/ p i n ] and the fraction due to double charges 
[/dbz], at 16% vacancy density. 
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FIG. 8: Application of the finite-size scaling of in-plane sus- 
ceptibility to estimate T c /J m 0.404 (common crossing point 
of the data) at 16% vacancy density at the vortex-on-vacancy 
critical anistropy. 
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FIG. 10: Out-of-plane susceptibilities \ zz vs - temperature 
for L = 128, at the three anisotropies studied. The lower 
curves at low T (open symbols) correspond to p vac = 0, the 
upper curves (solid symbols) correspond to p va c =0.16 . Un- 
like X° x or X vl \ there is only a very weak dependence of \ zz 
on L, mostly in the high-temperature phase. 



for the XY model. Comparing the results at the different 
anisotropies, there is no sudden change in the vacancy ef- 
fects, as far as can be seen from these limited data. The 
out-of-plane fluctuations vs. T for these nonzero A do not 
exhibit any particularly significant features due to the 
presence of vacancies. Generally, in the low-temperature 
phase, x zz increases with vacancy density, but even more 
so with increasing A, as summarized in Fig. [h]. It is clear 
that the out-of-plane fluctuations are aided by the pres- 
ence of vacancies, but from the limited data here, no 
significant conclusion about the role of pinned vortices 
vs. free vortices can be drawn. 



IV. DOUBLY-CHARGED VORTEX 
CONFIGURATIONS FROM SPIN RELAXATIONS 

Having seen the appearence, in general, of doubly 
charged vorticity localized on the vacancies, it is impor- 
tant to consider the basic analysis of their energetics. 
Clearly, in continuum theory, the static vortex energy 
(dependent on an integral of the form J J d 2 x |V0| 2 s» 
Jirq 2 lxi(R/a) ) is proportional to the squared charge. 
Therefore, one expects that the doubly-charged vortices, 
even when pinned on vacancies, should have considerably 
higher energy than singly-charged vortices (either pinned 
or free). Apparently, the absence of a spin at the center 
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FIG. 11: Various total system energies with a vortex present 
versus system radius R. 



of the q = 2 vortex, and the missing four interior bonds, 
significantly reduces the energy and allows for stability. 

Following the procedure in Ref. |?j, some doubly- 
charged vortex configurations were investigated numer- 
ically for their stability as a function of the anisotropy. 
For simplicity, a circular system of radius R, with sites 
defined on a square lattice, was used. The vacant site 
was placed exactly at the center of the circular system. 
Free boundary conditions applied at the edge of the sys- 
tem. The initial in-plane spin angles were set to that for 
a charge q vortex centered at position (x v ,y v ), 

<f ) ( x ,y) = q ttm- 1 (y^)+cf )0 . (16) 

For a q = ±2 vortex at the center of the system, a con- 
venient way to implement this expression on the xy spin 
components for arbitrary constant 4>o = 0, without using 
trigonometric functions is 

S*(x,y) = 4^4, Sy(x,y) = ^L, (17) 
x + y x + y 

the ± signs producing vortex/antivortex configurations. 

In order to test the in-plane to out-of-plane stability, 
all out-of-plane spin components were given small initial 
values S z — 0.001, thereby biasing the spin configura- 
tion possibly to go out-of-plane along the positive z-axis. 
After this small perturbation, all spins were normalized 
to unit length. The spin configuration was relaxed itera- 
tively by setting each spin to point along the direction of 
the effective field due to its neighbors, keeping the spin 
length fixed at unity. This leads eventually to a final con- 
figuration that is a local energy minimum of the Hamil- 
tonian, i.e., some form of stable configuration evolved 
from the original state, in some case with vorticity still 
present, and in other cases, no net vorticity. 



A. q = 2 vortex relaxation for XY model (A = 0) 

The first numerical relaxations were applied for the XY 
limit, A = 0, to get the general idea of the energy com- 
pared to that for q = 1 vortices. Typically, the energy 
found at A = should be expected to apply rather accu- 
rately to larger values of A, as long as the vortex remains 
in the planar configuration. These relaxations were per- 
formed for systems with radius ranging from R = 5a 
to R = 500a, as the energy is expected to have a loga- 
rithmic dependence on R. The energy results E vv for a 
q = 2 vortex centered on the vacancy are shown in Fig. 
|ll], and compared with similar results for q — 1 vortices. 
Additionally, the vortex energies E vp are shown when 
centered in a plaquette. For q = 1, this energy was found 
by relaxation to a stable vortex state, whereas, for q = 2, 
expression (|l^) was used to set the vortex centered in 
the plaquette, after which the energy was directly evalu- 
ated. This latter configuration for q = 2 is unstable, but 
was used for estimation of the vortex-on- vacancy binding 
energy, see below. 

Inspection of Fig. [ll] shows that, as expected, the 
doubly-charged vortices have considerably higher energy 
compared to singly charged, and furthermore, there is 
a nearly constant energy gap between the vortex-in- 
plaquette and vortex-on-vacancy states. Each data set 
fits extremely well to a logarithmic dependence on R in 
the form E = A+B \n(R/a). For q = 1, both curves have 
slope parameter B\ ss 3.17 J S 2 . For q = 2, both curves 
have slope parameter B2 ~ 12.7 J S 2 , a value very close 
to four times as large as that for q = 1, as might be ex- 
pected. The extra energy requirement for the q = 2 vor- 
tices clearly leads to a restriction on their thermal pop- 
ulation compared to q = 1 vortices. In all cases shown, 
the final spin configuration was found to be completely 
in-plane (all S z = 0). 

The difference between the vortex-on-vacancy and 
vortex-in-plaquette energies can be taken to define an 
energy for binding or pinning the vortex on the vacancy, 

AE q = E q ^ p — E q , vv . (18) 

Using the results shown, the binding energy for a q = 1 
vortex-on-vacancy is AEi ss 3.177JS' 2 , using the asymp- 
totic value as R — * 00. For doubly charged vortices, the 
binding energy is moderately higher, AE2 ~ 5.73JS" 2 , in 
contrast to the considerably higher creation energy for 
q = 2 vortices compared to q = 1 vortices. This result, 
however, must be taken with caution, since there is no 
actual stable q — 2 vortex free from a vacancy. 

An alternative view of the q = 2 vortex-on-vacancy 
energy would be to compare it to twice the energy of a 
system with a single q — 1 vortex centered in a plaque- 
tte, (2i5i vp ), because that is a stable state of the same 
total vorticity. However, the energy of the two q = 1 
vortices, in their own isolated systems, is always consid- 
erably less than that of a single q = 2 vortex. This is 
because 2Ei vp does not include the interaction potential 
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FIG. 12: After relaxation of a g = 2 vortex initially centered 
on an isolated vacancy in a circular system of radius R — 50, 
the total system energy (solid curve) is shown as a function 
of the anisotropy constant A. The vertical bars indicate the 
net out-of-plane magnetization of the relaxed configuration, 
on the same numerical scale. 



that would be present between two q = 1 vortices within 
the same system, which increases with the logarithm of 
their separation. Thus it is not a good reference number 
for estimation of the q — 2 binding energy on a vacancy. 



B. Anisotropy dependence of q = 2 vortex 
relaxation 



A preliminary analysis of the stability of a doubly- 
charged vortex can be performed by looking at the de- 
pendence of the relaxed configuration on the anisotropy 
parameter A > 0. It might be expected that a q = 2 
vortex could take on nonzero out-of-plane components 
when A becomes adequately close to 1, i.e., at weak easy- 
plane anisotropy, in a manner similar to the out-of-plane 
crossover for q = 1 vortices. The critical anisotropy could 
be expected to be different than the value A cv ~ 0.9545 
for pinned q = 1 vortices. Therefore, a limited num- 
ber of numerical experiments were realized on a circular 
system of radius R = 50a for various values of A above 
zero. Again, the initial condition was ag=2 vortex cen- 
tered on the vacancy at the center of the system, with 
small positive out-of-plane components (S z = 0.001) at 
all sites. 

Certain aspects of these results are summarized in Fig. 
|l2| , where the energy of the state obtained after the relax- 
ation is plotted versus the anisotropy parameter A that 
was used. There are several types of results, depending 
on the range of A being considered. 

For the whole range < A < 0.545 (region P2), the 
isolated q = 2 vortex remains in a stable planar con- 
figuration on the vacancy, with no out-of-plane mag- 



netization, and relatively high energy. For the narrow 
range 0.545 < A << 0.57 (region 2 ), the q — 2 vor- 
tex remains stable on the vacancy, but develops nonzero 
out-of-plane magnetization, with an insignificant reduc- 
tion in energy. The net out-of-plane magnetization of 
the relaxed state is indicated in Fig. 12 by the bars ex- 
tending above the energy curve. Total M z grows until 
A reaches about 0.57, at which point the q — 2 vortic- 
ity concentrated on the vacancy becomes unstable, and 
breaks into one q — 1 in-plane vortex on the vacancy, 
and a nearby free 5 = 1 in-plane vortex. This situation 
holds for 0.57 < A < 0.66 (region PP); the configura- 
tion has zero out-of-plane magnetization once again, and 
lower energy than that for the q = 2 vortex pinned on a 
vacancy. As A increases within this range, the free vortex 
progressively moves farther from the vacancy. 

When A ranges from about 0.67 to 0.68 (region PO), 
the free vortex starts to develop a nonzero positive out- 
of-plane component, while the pinned vortex remains pla- 
nar. Finally, at A w 0.685 and above (region OO), the 
relaxed configuration consists of two q — 1 positively po- 
larized out-of-plane vortices centered symmetrically on 
opposite sides of the vacancy. For example, the relaxed 
configuration obtained for A = 0.7 is shown in Fig. [b]. 
As A increases, the separation of the pair increases at the 
same time that their out-of-plane component increases, 
while the energy decreases. Eventually, the separation 
surpasses the diameter of the system, and the vorticity 
escapes out the boundary, leaving a final configuration 
with uniform magnetization and zero energy. This oc- 
cured for A > 0.77 in the system of radius R = 50. 

It is apparent that a localized q = 2 vorticity has 
severely limited stability, compared to q = 1 vortices. 
The q = 2 vorticity even tends to grow out-of-plane com- 
ponents as a way to enhance its stability, but this has a 
very limited range of utility (region C^)- Once the vortic- 
ity splits into individual q = 1 vortices, they are seen to 
influence each other, probably via an interaction with the 
vacancy. This is apparent because out-of-plane compo- 
nents begin forming for A below the critical anisotropy 
parameter A c relevant for vortices far from vacancies. 
Furthermore, the pairs of out-of-plane q — 1 vortices in 
region OO appear to repel each other, while at the same 
time being attracted to the vacancy, which would lead 
to a mechanically stabilized configuration. Inspection of 
the spin configurations in region OO (as in Fig. |l3|) shows 
spin components of one vortex to be completely symmet- 
rical to the spin components in the other vortex, when 
reflected across the center of the system. 



V. CONCLUSIONS 

In the Monte Carlo and spin dynamics calculations 
presented here for a 2D easy-plane anisotropic Heisen- 
berg model, the presence of vacancies has been seen to 
affect the details of the BKT transition and the types of 
vorticity present. 
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FIG. 13: Final state of relaxation of a q — 2 vortex initially 
centered on an isolated vacancy in a circular system of radius 
R — 50 with A = 0.7 (only the central region of system is 
shown). Part (a) shows the projection of xy spin components 
on the plane, as explained in Fig. ^j. In part (b) the lengths 
of the arrows are equal to the 2 spin components, while the 
directions are still given by the xy projections. 



As seen in the planar rotator model, T c is lowered by 
the presence of vacancies. This naturally results because 
the disorder in the transition becomes dominated by the 
generation of vortices pinned on the vacancies. When 
formed centered on vacancies, the vortex energy is signifi- 
cantly lower than that for vortices centered in plaquettes. 



FIG. 14: A spin configuration for L = 32, A = 0, p vac = 0.16, 
at T ~ T c ~ 0.48J, where 9 out of the 10 vortices present 
have formed on vacancies. 



Indeed, 9=1 vortices pinned on vacancies in a square lat- 
tice, have a formation energy of about 3.17 J S 2 less than 
when centered in plaquettes, while the transition tem- 
perature corresponds to an energy less than US 2 . Thus, 
at temperatures near and below T c , the small amount 
of vorticity that is present is predominantly pinned on 
vacancies, such as in Fig. [l4| These pinned vortices are 
initiating and controlling the transition. The vacancies 
are the nucleation sites for the spin disordering. On the 
other hand, vacancies reduce the rate at which total vor- 
ticity density rises in the high-temperature phase (Figs. 

ID- 

At larger anisotropy parameter A c , it is known that the 
vortex-on-vacancw energy is much closer to the vortex-in- 
plaquette energy.Q For example, at A = 0.99, the differ- 
ence in these energies is only 0.23JS 2 . Then one might 
expect a lesser dominance of pinned vorticity, however, 
that does not appear to be the case at A = A cv - There is 
no substantial qualitative change in the fraction of pinned 
vorticity when compared to the XY model. Qualitatively 
speaking, the details of the BKT transition at higher A, 
with vacancies, are not significantly different than those 
found for the XY model. 

The presence of vacancies leads to a new effect, namely, 
the generation of doubly-charged vorticity that is sta- 
ble when centered on vacancies. In thermal equilibrium, 
this effect apparently occurs regardless of the easy-plane 
anisotropy strength. In general, these would be thermo- 
dynamically prohibited, due to their larger energy, based 
on the usual dependence of vortex creation energy on 
charge squared. They still have significantly higher en- 
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ergy than singly charged vortices centered on vacancies, 
however, the missing bonds in the core region help to 
reduce their total energy compared to what they would 
have if centered in a plaquette. 

A spin-relaxation energy minimization shows that an 
individual doubly-charged vortex centered on a vacancy 
in a circular system may be stable only for a limited range 
of anisotropy constant. The q = 2 vortex-on- vacancy 
stays stable and planar for < A < 0.545. In a very 
narrow range, 0.545 < A < 0.57, the q = 2 vortex-on- 
vacancy still remains stable, but with a small out-of-plane 
component. For A > 0.57, it does not appear to be stable, 
but instead breaks apart into two q = 1 vortices that re- 
pel each other while being attracted to the vacancy. One 
might define a lower critical anisotropy A cv ,i ~ 0.545 for 
the in-plane to out-of-plane q = 2 stability, and an upper 
critical anisotropy A cv ,2 ~ 0.57 for the breakdown into 
lower charged vorticies. In contrast, there is no choice of 
anisotropy constant that stabilizes a, q = 2 vortex in the 
center of a plaquette. 

These results are intriguing, because even though they 
show a limited range of stability for doubly-charged vor- 
ticity, nevertheless, these excitations appear in the MC 



simulations at A = 0.7 and A = 0.9545, above the criti- 
cal anisotropy parameters. Of course, one could always 
search groups of four plaquettes in the pure model also 
to find localized vorticity of double charge (two q = 1 
vortices in neighboring plaquettes), although it would 
appear very rarely, due to the mutual repulsion of the 
vortices. The difference here, is that the presence of 
a vacancy attracts vorticity and certainly enhances the 
chances to find doubled vorticity within the area of four 
neighboring plaquettes. In addition, the spin relaxations 
show that the doubly charged vortex can be a static ob- 
ject, which can never be expected for q = 1 vortices in 
neighboring fully occupied plaquettes. 
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